Solitons, solitonic vortices, and vortex rings in a cylindrical Bose-Einstein condensate 
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Quasi-one-dimensional solitons that may occur in an elongated Bose-Einstein condensate become 
unstable at high particle density. We study two basic modes of instability and the corresponding 
bifurcations to genuinely three-dimensional solitary waves such as axisymmetric vortex rings and 
non-axisymmetric solitonic vortices. We calculate the profiles of the above structures and examine 
their dependence on the velocity of propagation along a cylindrical trap. At sufficiently high ve- 
locity, both the vortex ring and the solitonic vortex transform into an axisymmetric soliton. We 
also calculate the energy-momentum dispersions and show that a Lieb-type mode appears in the 
excitation spectrum for all particle densities. 

PACS numbers: 03.75.Lm, 05.30.Jp, 05.45.Yv 



I. INTRODUCTION 

The possible occurrence of solitons in a Bose-Einstein 
condensate (BEC) was theoretically predicted sometime 
ago, within a one-dimensional (ID) Gross-Pitaevskii 
model P, , but experimental observation was delayed 
because of the absence of a physical realization of a 
strictly ID Bose gas. Nonetheless, similar coherent struc- 
tures were recently observed in BECs confined in traps 
of varying geometry 0, 3 by a method (phase imprint- 
ing) that was directly inspired by the analytical structure 
of ID solitons. Yet it has become clear that quasi-lD 
solitons are susceptible to various instabilities within the 
three-dimensional (3D) environment of realistic traps. 

Indeed, a stability analysis based on the linear 
Bogoliubov-de Gennes (BdG) equations jsf revealed that 
an axisymmetric soliton is stable only at sufficiently low 
particle number or high aspect ratio in an elongated trap. 
The primary mode of instability was shown to result 
from non-axisymmetric perturbations, in the sense that a 
purely imaginary eigenvalue first appears in the spectrum 
of the BdG equations with azimuthal angular momentum 
m = l. This "snake instability" was further analyzed in 
Ref. and argued to be responsible for a possible decay 
of a soliton into vortex rings and/or vortices. In fact, 
based on the above theoretical work, non-stationary vor- 
tex rings were experimentally detected in a spherical trap 
0. On the other hand. Brand and Reinhardt [1,0 sug- 
gested that the primitive mode associated with the snake 
instability is a stable structure that may be called soli- 
tonic vortex. 

The emerging picture is sufficiently perple xing to de- 
serve closer attention. In some recent work [l^ [ll| we 
have carried out a detailed calculation of axisymmetric 
solitons in a cylindrical trap, which were shown to be- 
come unstable at high particle density even if we restrict 
ourselves to radial (m = 0) perturbations. The soliton 
was then shown to bifurcate to an axisymmetric vortex 
ring with lower energy. Here, we complete this work by 
allowing azimuthal deformations of the soliton, in order 



to eventually account for the primary (m = 1) instabil- 
ity discussed in Ref. 0, Q ■ We thus obtain a reasonably 
complete description of the two basic modes of instability 
of the soliton and the corresponding bifurcations to ax- 
isymmetric vortex rings and non-axisymmetric solitonic 
vortices. 

The simplest theoretical picture is obtained in the ideal 
limit of an infinitely elongated cylindrical trap, because 
static solitary waves may then be treated on equal foot- 
ing with those propagating at nonzero velocity. Hence, 
all calculations presented in the main text pertain to 
an infinite cylindrical trap. In Sec. II we formulate the 
problem and briefly describe the numerical methods em- 
ployed. The three basic types of static solitary waves 
(solitons, solitonic vortices, and vortex rings) and their 
relative stability are discussed in Sec. Ill, whereas soli- 
tary waves moving along the cylindrical trap at constant 
velocity are calculated in Sec. IV. The main conclusions 
are summarized in Sec. V, and some important technical 
details are relegated to two appendices. In particular, 
the case of a finite axisymmetric trap is discussed in Ap- 
pendix B. 



II. FORMULATION 

An axisymmetric harmonic trap is characterized by a 
transverse confinement frequency lo^_ and a longitudinal 
frequency wy. We first consider the special limit of a 
cylindrical trap, with uj\\ = 0, a restriction that is re- 
laxed only in Appendix B. In this limit, physical units 
are rationalized as in Ref. pTj. Thus time is measured 
in units of l/uj±_ and distance in units of the transverse 
oscillator length = (Ti/Mw^)-'^/^ where M is the mass 
of each atom. Complete specification of the system also 
requires as input the linear density u which is equal to 
the average number of particles per unit length of the 
cylindrical trap. An important dimensionless parameter 
is then given by ^ = va where a is the scattering length. 
A field rescaled according to — > ^Jv^ ja^^ satisfies the 
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rationalized Gross-Pitaevskii equation 
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(1) 



where T4r = {x'^ + J/^)/2 = is the trapping potential 
and g = Att"/ is a dimcnsionless coupling constant that is 
the only free parameter. The conserved energy functional 
associated with Eq. is given by 



W 



\V^\^ + 2VtA^\^+g\^\^)dV (2) 



and yields energy in units of 'yj_(hu}±) where 7^ = va± = 
l{a±/a). 

We have also considered various forms of dissipative 
dynamics, mainly as relaxation algorithms for the calcu- 
lation of stationary states. The simplest possibility is to 
make the replacement i d'^/dt {i—Qd'^/dt in Eq. J^l, 
where C, is some dissipation constant. An extreme special 
case widely employed in numerical simulations is the 
fuUy-dissipative dynamics obtained by the replacement 
id'^/dt — > —d'^/dt. However, in either case, particle 
number conservation is violated and special precautions 
are necessary to enforce a definite number of particles. 

An interesting alternative that preserves particle num- 
ber was originally introduced by Carlson for the 
study of the homogeneous BEC and is simply adapted 
to the case of a confined BEC. Actually, we have em- 
ployed a special case of Carlson's dissipative dynamics 
governed by 
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dt 



A* + Vtr* + <7l*P*- (3) 



An immediate consequence of Eq. iPJ is the continuity 
equation 



dn 
'dt 



V- J = 0. 



(4) 



irrespective of the specific value of the dissipation con- 
stant C- Here, ri=|^'p is the local particle density and 



(5) 



is the usual current density, where </> is the phase of the 
wave function \1/ = -^71 e"^ . Therefore, the total number 
of particles is indeed conserved by the evolution equa- 
tion |j2Jl. However, the total energy defined from Eq. Q 
dissipates according to 



dW_ r /dn 
'dr~~. [~dt 



dV, 



(6) 



cubic lattice to test the stability of solitons and even- 
tually obtain rough estimates of other possible station- 
ary states. These estimates are then fed into a Newton- 
Raphson algorithm of the type considered in Ref. to 
obtain our final results more accurately as well as effi- 
ciently. The latter algorithm is currently upgraded to a 
fully-3D iterative scheme that does not necessarily en- 
force axial symmetry. Explicit results are discussed in 
the following sections. 



III. STATIC SOLITARY WAVES 

We first consider a special class of stationary states 
obtained by making the substitution 'I'(r, t) e^'^*^'(r) 
in Eq. to derive the static equation 



-A* 



1 



(7) 



where fi is the chemical potential (in units of huj±) re- 
quired to enforce a definite number of atoms. The sim- 
plest and most important solution of Eq. Q is the ground 
state wave function ^' = ^'o (p) which depends only on the 
radial distance p= (x^ -I- j/^)^^^ and may be chosen to be 
real and positive. It also satisfies the constraint 



27: pdp-^lip) = 1, 



(8) 



by a specific choice of the chemical potential p, = p{j) 
for each value of the dimcnsionless coupling constant 7. 
Explicit illustrations of the ground state may be found 
in Ref. (111] . 

A more general class of excited self-consistent states is 
obtained by solving Eq. (|7|) with the same chemical po- 
tential, as is appropriate on an infinitely elongated trap, 
and boundary conditions 



lim \'i>(x,y,z)\ 

-^±00 



9* 



0. 



(9) 



The class of solutions actually discussed in the present 
paper is further restricted by the parity relation 



'^{x,y, -z) = **(a;,y,z) 



(10) 



modulo an overall constant phase that is set equal to 
zero. Relation UlUI) is obviously consistent with Eq. iQ. 

To conclude the general description, we also consider 
the excitation energy defined as 

E = W -W^-p [ - l^-oH dV , (11) 



as long as the particle density n remains time dependent. 

Although both types of dissipation discussed above are 
employed in our numerical calculations, Carlson's pro- 
posal has some advantages for our current purposes. We 
have thus developed a discrete version of Eq. (3) on a 



where both W and Wq are calculated from the energy 
functional ^ applied for a general self-consistent state 
and the ground state ^'o, respectively. The presence 
in Eq. (|11|) of a term that is proportional to the chemical 
potential p ensures that we are calculating the energy 



FIG. 2: Profile of a static {v = 0) solitonic vortex for 7 = 
3, illustrated through density plots over two planes that are 
perpendicular to each other and both contain the symmetry 
(2) axis of the cylindrical trap. Regions with high particle 

^ ,1,1,1,1, density are bright while regions with zero density are black. 

2 4 6 8 10 The size of each frame is [-5,5] X [-5,5]. 
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FIG. 1; Excitation energy E in units of 'y_L{hu)_i_) as a function 
of the dimensionless coupling constant 7, for static solitary 
waves such as a soliton (S), a solitonic vortex (SV), and a 
vortex ring (VR). Bifurcations occur at the two critical cou- 
plings 70 = 1.5 and 71 = 4. 

difference between two states with the same number of 
particles. 

In recent work 0, we restricted ourselves to ax- 
isymmetric states described by a wave function of the 
form 'i' = '^{p,z). The simplest possibility is a static 
(black) soUton whose wave function is odd under the par- 
ity reflection z ^ —z, and purely imaginary in view of 
the phase convention adopted in Eq. H1U|I . The excita- 
tion energy of a static soliton calculated from Eq. 1)11(1 
is illustrated in Fig. ^ as a function of the dimensionless 
coupling constant 7. 

When 7 exceeds the critical value 71 =4, a new class of 
axisymmetric solitary waves emerges with the standard 
characteristics of vortex rings |3 • Unlike the case 
of a homogeneous BEC, a vortex ring may now be static 
and is described by a complex wave function that satisfies 
the parity relation ^'(p, —z) — z). The calculated 

excitation energy of a static vortex ring is also shown in 
Fig. n and is seen to smoothly bifurcate from the soliton 
branch at the critical coupling 71 , while it remains consis- 
tently lower than the soliton energy for 7 > 71 . Therefore, 
one may conclude that a radial (m = 0) perturbation of a 
soliton will render it unstable when 7 > 71 and eventually 
transform it into an axisymmetric vortex ring. 

One is tempted to presume that the soliton is stable for 
7 < 71. However, the linear stability analysis of Ref. 0,0 
indicates that the primary soliton instability sets in at a 
lower critical coupling 7o<7i through azimuthal (m=l) 
perturbations. To probe for such a possibility, we solve 
the dissipative Eq. Q with initial condition supplied by 
the calculated static soliton which is found to remain sta- 



ble at all times when 7 < 70 = 1.5. However, for 7 > 70, the 
inherent numerical noise introduces an instability that 
transforms the soliton into a non-axisymmetric structure 
which is subsequently rectified by feeding it into a 3D 
Newton-Raphson iterative algorithm that is designed to 
directly solve the static equation Q. 

We thus obtain a static solitary wave that is partially 
illustrated in Fig.|21for 7 = 3 through density plots over 
two planes that are perpendicular to each other and both 
contain the z axis. We have furthermore examined the 
phase to establish that the calculated wave function de- 
scribes a deformed vortex that pierces through the trap 
in a direction that is perpendicular to its symmetry (z) 
axis. In Fig. [3 the vortex axis is identified with the x 
axis as a matter of convention, whereas a degenerate class 
of equivalent configurations may be obtained by simple 
azimuthal rotations around the z axis. The resulting soli- 
tary wave displays the characteristics of a solitonic vortex 
calculated earlier by Brand and Reinhardt 0, Q in traps 
of varying geometry. 

The energy of the solitonic vortex is depicted in Fig.Q] 
as a function of 7 and is seen to emerge from the soli- 
ton branch at the critical coupling 70 = 1.5, while it re- 
mains smaller than the energy of both the soliton and the 
vortex ring for all 7 > 70- The physical picture can be 
explained by returning to the initial-value problem asso- 
ciated with the dissipative Eq. ||3Jl using the static soliton 
as initial condition. As mentioned already, the soliton re- 
mains stable for 7 < 70 but decays into a solitonic vortex 
when 7o < 7 < 71 . The pattern is more complex for 7 > 71 
where the soliton initially deforms radially to become an 
axisymmetric vortex ring which eventually decays into a 
stable solitonic vortex. 

Therefore, a vortex ring is strictly speaking unstable. 
But our real-time simulations indicate that vortex rings 
are sufficiently long lived to be observed, as was actu- 
ally the case in the experiment of Ref. ^ - albeit in the 



special limit of a spherical trap. On the other hand, it 
is curious that the predicted solitonic vortex has not yet 
been seen in experiment. 

The picture is completed with a few remarks. One 
should expect that there exists a series of additional crit- 
ical couplings (densities) at which the soliton bifurcates 
into multiple vortex rings and/or solitonic vortices. For 
example, a double vortex ring was calculated in Ref. [Tlj 
for 7 > 72 = 12. In other words, we have only described 
the two basic patterns of instability that correspond to 
the appearance of unstable modes in the m = and 
m = 1 sectors of the linear BdG equations. These two 
modes may be dominant at intermediate densities, but 
more complicated patterns should be expected to occur 
at higher densities, as is evident in the experiment of 
Ref. 7] which is further discussed in Appendix B. 



IV. MOVING SOLITARY WAVES 

All three types of static solitary waves discussed in the 
preceding section are special members of a more general 
class of stationary states that propagate along the cylin- 
drical trap with constant velocity v. These are obtained 
by making the substitution ^'(r,i) e~*''*^'(a;, y, ^) in 
Eq. 1^ , with — z — vt, to derive the stationary differen- 
tial equation 



V = 0.4 



+ lp^^5> + g\^\H , (12) 



A = 



de 



which reduces to the static equation {Tj) at w = 0. The 
boundary conditions ® and the parity relation (fTU]) are 
also generalized by the simple replacement z — > ^. 

Equation l|12|) is solved by feeding a static {v = 0) soli- 
tary wave into the Newton-Raphson algorithm and by 
upgrading the velocity typically in steps of 5v = ±0.1. 
The results for axisymmetric solitons and vortex rings 
already presented in Ref. 10., 11] have been confirmed 
using the currently employed fully-3D algorithm. There- 
fore, in the following, we shall mostly discuss the case 
of a moving solitonic vortex for which no calculation has 
been given in the literature. The static solitonic vortex 
of Fig. 121 is significantly reshaped when it moves, as is 
illustrated in Fig.|3|for two values of the velocity (?; = 0.4 
and 0.8) that are both smaller than the speed of sound 
c= 1.3 calculated for 7 = 3 [ij. The w = 0.4 entry of Fig.EI 
makes it clear that a moving solitonic vortex is shifted off 
center. While the vortex axis remains within the basal xy 
plane, it is clearly displaced away from the center of the 
trap. The displacement increases with increasing speed 
until we reach a critical speed vq = 0.8 where the vortex 
disappears from the trap, in the sense that the solitary 
wave is converted into a gray axisymmetric soliton. This 
implies, in particular, that the density plots over the xz 




V = 0.8 




y 



FIG. 3: Profile of a moving solitonic vortex for 7 = 3 and two 
values of the velocity v = 0.4 and 0.8. Conventions are the 
same as in Fig. |21 except that ^ = z — vt now measures the 
distance from the center of the moving solitonic vortex. The 
speed of sound calculated at 7 = 8 is c=1.3. 



and yz planes become identical, as is apparent in the 
v = 0.8 entry of Fig. El 

For vo <\v\ <c the calculated solitary wave is indeed 
identical to a moving axisymmetric soliton and reduces to 
a weakly nonlinear soundlike pulse in the limit \v\ — > c. 
A corollary of this discussion is that a high-speed soli- 
ton can be stable even for couplings 7 that are greater 
than the critical coupling 70 = 1.5, in apparent agreement 
with a similar conclusion reached in Ref. that is now 
made more precise. Specifically, for 7 < 70 a solitonic vor- 
tex does not exist and the soliton is stable over the entire 
range of velocities < |w| < c, as expected. For 7 > 70, 
the soliton is stable over the limited range vq < \v\ < c, 
where vq = ^0(7) is the critical speed at which the soli- 
tonic vortex is converted into a soliton and = 0(7) is the 
calculated speed of sound Clearly, tig(7 = 79) = 0, 

while a numerical investigation of the ratio vo{'j)/c{'j) 
shows that it approaches unity with increasing 7. How- 
ever, we are not in a position to ascertain whether or 
not there exists a finite (critical) 7 above which vq = c 
and, hence, the soliton is not stable for any velocity. In 
any case, the solitonic vortex is the fundamental mode 
as long as 7 > 70 . 

Further insight into the nature of a moving solitonic 
vortex is obtained by considering the behavior of its 
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FIG. 4: Asymptotic phase difference as a function of velocity 
for a solitonic vortex (SV, solid line) and a soliton (S, dashed 
line) . The two curves merge at the critical velocities v — itvo 
where uo = 0.6 c. 



phase (j) — (/>(a;,y,^) in the asymptotic limits ^ — > ±oo. 
Let (p± — 4>{x,y,£, —±00) be the two asymptotic values 
of the phase, which are, in principle, some functions of 
X and y. However, detailed numerical evidence suggests 
that 0+ and (j)^ are, in fact, constants that depend only 
on the velocity v. This conclusion may be substantiated 
in part by an analytical argument briefly described in 
Appendix A. Therefore, the asymptotic phase difference 

,50 = 0+ - 0_ (13) 

depends only on the velocity v and characterizes a soli- 
tary wave. Incidentally, the above conclusion applies to 
all three types of solitary waves discussed in the present 
paper and is pertinent to their experimental production 
through phase imprinting 0, 0| . 

As an illustration, we return to the solitonic vortex 
calculated for 7 = 3 and plot —5(j) as a function of velocity 
by a solid line in Fig.^ whereas the corresponding result 
for an axisymmetric soliton at 7 = 3 is shown by a dashed 
line. In both cases, —S(j) = 7T for v = 0, and —(50 = (mod 
27r) in the extreme limits v — > ±c which are not actually 
reached in Fig. 0] due to numerical inaccuracies. The 
same figure provides clear evidence for the conversion of 
a solitonic vortex into an axisymmetric soliton at the 
critical speed wo = 0.8 or wo/c = 0.6. 

The preceding discussion of the asymptotic phase is 
also pertinent to the definition of the impulse of a solitary 
wave pH which is now simplified to 

Q = V-S(j), V=fn^dV, (14) 
J oz 

where V is the usual linear momentum and (50 is given 




FIG. 5: Excitation energy E in units of 7x(?i(jJ±) vs impulse 
Q in units of ^±{h/a±) for a soliton (S) at 7 = 1. 



by Eq. (|13|l . The impulse Q is measured in units of hv = 
j±{h/a±). 

A solitary wave may then be characterized by the im- 
pulse Q and the excitation energy E defined in Eq. Hll(l . 
Both E and Q are definite functions of the velocity v that 
can be eliminated to yield the dispersion E = E{Q). A 
nontrivial check of consistency is provided by the group 
velocity relation v — dE/dQ which was accurately con- 
firmed by our numerical calculation. For comparison 
purposes, we reproduce in Fig. [S] the dispersion of a sta- 
ble axisy mmetric soliton calculated for 7 = 1 < 70 in 
Ref. 11]. This dispersion displays the characteristics of 
a mode originally derived by Lieb 14] within a full quan- 
tum treatment of a homogeneous ID Bose gas based on 
the Bethe Ansatz, and was later rederived by a semiclas- 
sical calculation based on the solitary wave of the ID 
Gross-Pitaevskii model 0, . A similar mode was also 
derived within an effective ID model that is designed to 
approximately describe a cylindrical trap 0| . 

We are now in a position to calculate the dispersion 
of a moving solitonic vortex. The dispersion is shown 
by a solid line in Fig. |H1 and is seen to retain the basic 
characteristics of a Lieb-type mode. We further plot by 
a dashed line the dispersion of an axisymmetric soliton 
which is expected to be unstable at 7 = 3. However, 
the two dispersions shown in Fig. El merge at two critical 
points that correspond to velocities v = ±vq where vq = 
0.6 c is the critical speed at which the solitonic vortex is 
converted into an axisymmetric soliton. 

A subtle detail not seen in our limited illustrations is 
that the calculated soliton for 7 = 3 exhibits some ringlike 
features, when v ^ 0, which are not present for v = 
and also disappear before reaching the critical speed i^o- 
These features could be viewed as preliminary signals for 
the formation of fully-fledged vortex rings that occur for 
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FIG. 6: Excitation energy E vs impulse Q for a solitonic 
vortex (SV, solid line) and a soliton (S, dashed line) at 7 = 3. 
The two dispersions merge at two points that correspond to 
velocities v — iizvo where no = 0.6 c. 



FIG. 7: Excitation energy E vs impulse Q for a solitonic 
vortex (SV, solid line), a vortex ring (VR, dotted dashed line) 
and a soliton (S, dashed line) at 7 = 7. The SV and VR 
dispersions merge at two points that correspond to velocities 
v = ±vo where wo=0.8c. 



stronger couplings in the region 7>7i=4. 

Since a detailed calculation of both static and moving 
vortex rings was already described in Ref. 0, we 
merely complete the picture by illustrating in Fig. [Tithe 
dispersions of a soliton, a vortex ring, as well as a soli- 
tonic vortex calculated for 7 = 7 where all three possibil- 
ities occur according to the bifurcation pattern of Fig.Q 
The solitonic vortex is still characterized by a Lieb-type 
dispersion that merges with the dispersion of the vortex 
ring at the critical speed uo = 0.8c, where c=1.6 is the 
speed of sound calculated at 7 = 7. For vo<\v\<c both 
the solitonic vortex and the vortex ring become identical 
to a gray axisymmetric soliton. In fact, the vortex ring 
transforms into a soliton at a speed somewhat smaller 
than vq. 

Therefore, it is fair to say that the primary mode as- 
sociated with the snake instability of a soliton is indeed 
the solitonic vortex Q whose calculation has been sig- 
nificantly extended in the present paper to account for 
steady motion along a cylindrical trap. Actually, one 
may envisage a further generalization to a spinning soli- 
tonic vortex that undergoes rigid-body precession around 
the z axis at constant frequency u>, while it moves along 
the same axis at constant velocity v. One could then 
calculate the excitation energy E = E(lu,v), the im- 
pulse Q — Q(u;,v), and the azimuthal angular momen- 
tum Lz — Lz(u},v) to eventually derive a dispersion of 
the form E — E{Q,Lz). The Lieb-type dispersion of a 
solitonic vortex shown by a solid line in Figs. and [7| 
would then be the lower edge of a continuum of states 
parametrized by the angular momentum L^- This in- 
teresting possibility will not be discussed further in the 



present paper. 



V. CONCLUSION 

A phase imprint Scj) applied between the two ends of an 
elongated trap leads to the formation of pure quasi- ID 
solitons only for weak couplings (densities) in the region 
7 < 7o = 1.5. For intermediate couplings in the region 
7o < 7 < 71 = 4, solitons may be produced temporar- 
ily but they eventually decay into stable solitonic vor- 
tices. A new pattern emerges for 7 > 71 where a phase 
imprint should initially produce unstable axisymmetric 
solitons that subsequently deform radially to become ax- 
isymmetric vortex rings and finally decay into solitonic 
vortices. This pattern is further complicated by the fact 
that high-speed vortex rings or solitonic vortices are in- 
distinguishable from axisymmetric solitons. 

It is thus important to reexamine the parameters of 
the experiment described in Ref. According to the 
estimates obtained in Appendix B, the length of the trap 
is significantly larger than its radius. Hence the approxi- 
mation by an infinite cylindrical trap may be reasonable. 
Furthermore, a conservative estimate of the effective cou- 
pling leads to 7 ~ 7 which lies in the region 7 > 71. 
Therefore, all three types of solitary waves theoretically 
discussed in the present paper should have been produced 
in the original experiment of Ref. 0. Also note that 
non-stationary vortex rings were actually detected in the 
experiment of Ref. but on a spherical trap that is sig- 
nificantly different than the cylindrical trap considered 
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here. Nevertheless, our present analysis can be extended 
to the case of a finite trap, as discussed in Appendix B. 

Finally, the calculation of the present paper clearly 
demonstrates that a Lieb-type mode appears in the ex- 
citation spectrum for all densities, either in the form of 
a quasi-lD soliton (for 7 < 70) or in the form of a soli- 
tonic vortex (for 7>7o). It is thus interesting to examine 
whether or not such a mode can be measured. 



Acknowledgments 

We are grateful to J. Brand for a stimulating discus- 
sion that brought the work of Ref. Q to our attention, 
to N.R. Cooper for his constant interest in the present 
effort, and to A.S. Fokas, and N.S. Manton for a num- 
ber of useful comments. This work was supported by the 
EPSRC under contract GR/R96026/01 (S.K.). 



APPENDIX A: ASYMPTOTIC PHASE 

The essence of the analytical argument given below 
was suggested to us by Fokas For a solitary wave 

traveling with constant velocity v, the continuity equa- 
tion ^ may be written as —vdn/d^ + V • (nV(/)) = 0. In 
the hmits ^— s-ioo, no{p) — \'^o{p)\'^ and all ^ (or z) 
derivatives vanish on account of the boundary conditions 
©. Hence, Vj_ • [no{p)^ j_(l>±{x, y)]=0, where Vj^ is the 
2D gradient operator and 4>±{x,y)=^(l){x,y,£,—±(X)) are 
the unknown asymptotic limits of the phase. If we further 
employ cylindrical coordinates {x — p cost/?, y — p soup) the 
asymptotic form of the continuity equation reads 



d 
dp 



pno{p) ■ 



dp 



m{p) 



(Al) 



For axisymmetric solitons and vortex rings — (j3± (p) 
and thus the second term in Eq. (|A1I) vanishes. A first 
integral is then given by d4>±/dp = ci/pno{p) where ci 
is some integration constant. Also note that the energy 
functional Q contains a term n{d(l)/dp)^ that becomes 
no{p){d(l)±/dp)^ = cl/p'^no{p) at the two ends of the trap. 
Since the ground-state density no{p) is finite at p = and 
vanishes in the limit p 00, the above term would lead 
to a badly divergent integral over the radial distance p 
unless the integration constant ci vanish. Hence, for an 
axisymmetric solitary wave with finite excitation energy, 
d(j)± /dp — and thus the asymptotic phases (f)± can de- 
pend only on the velocity v. 

On the other hand, we have not been able to extend the 
above argument to the case of non-axisymmetric solitary 
waves for which both terms in Eq. ljAl|l are potentially 
relevant. Nevertheless, our numerical results for the soli- 
tonic vortex discussed in the main text again suggest that 
the asymptotic phases 0± = (t>±ix,y) are, in fact, inde- 
pendent of x and y. 



APPENDIX B: FINITE TRAPS 

A finite harmonic trap with axial symmetry is char- 
acterized by the two confinement frequencies and 
cj|| . We now introduce rationalized quantities that are 
in some respects different than those employed for the 
discussion of the infinite cylindrical trap. Time and dis- 
tance are still measured in units of 1 /uj± and a± but the 
wave function is rescaled according to 5* ^ ViV^'/aj^^j 
where N is the total number of atoms, and thus satis- 
fies the constraint / j^'pdF — 1. The Gross-Pitacvskii 
equation is then formally identical to Eq. but 



47ra, Vtr 



1 



ip' + P'z'), 



(Bl) 



where a — Na/aj_ and f3 = uji^^/uj^ are now the two in- 
dependent dimensionless parameters. To make contact 
with the infinite trap, we calculate the average linear 
density from i'{z) — {N/a±) / j^'pdxdy and define the 
dimensionless quantity 7(2:) — ai'{z) or 



7(z) 



I^Pf dxdy , 



(B2) 



which is the analog of the dimensionless coupling con- 
stant 7 employed in the main text, except that it now 
depends on z. Nevertheless, 7(2:) becomes increasingly 
uniform in the limit a ^ 00 and /? — *■ 0, holding 7(0) = 7 
fixed, which is the ideal limit of an infinite cylindrical trap 
studied in the main body of the paper. 

As a simple illustration, we consider the Thomas-Fermi 
(TF) approximation of the ground state to find that 



7(2:) 



1 

16 



(B3) 



where the chemical potential is given by 2p— {15a(3)^/'^. 
Therefore, the estimated maximum value of 7 is given by 
7max = 7(0) = A more conservative estimate is the 
average value defined from ^f^^ — ajL where L = 2^^/2p/ (3 
is the TF length of the trap and thus 7av = (8/15)7niax- 
The TF radius of the trap is accordingly given by i? = 
\/2p,. For the parameters of the experiment described in 
Ref. j3l we find i? = 2 /z m and L = 115 m, while 7max = 
13 and 7av = 7. These estimates are further discussed in 
the concluding section in the main text. 

Now, static solitons, solitonic vortices, and vortex rings 
exist also in a finite trap and we have calculated them 
explicitly for various values of a and (3 to confirm the 
general picture described in Sec. III. Here, we briefly 
summarize the basic facts. For each value of the aspect 
ratio /3, there exist two critical couplings ao = ao(/3) and 
OLi = ai {P) that lead to a bifurcation pattern similar to 
that shown in Fig. ^ The critical coupling ao coincides 
with the coupling where a purely imaginary eigenvalue 
first appears in the m = 1 sector of the linear BdG equa- 
tions 5]. Similarly, the critical coupling ai is related to 
an instability that occurs in the m = sector, but the 
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FIG. 8: Profiles of a static soliton (left panel) and a static 
vortex ring (right panel) on a finite trap with a = 100 and 
/3 = 1/4, illustrated through density plots over a plane that 
contains the symmetry {z) axis and cuts across the axisym- 
metric trap. The complete pictures can be obtained by simple 
revolution around the z axis because both the soliton and the 
vortex ring are axisymmetric. 

picture is now significantly complicated by the appear- 
ance of complex modes . Our numerical solution of the 
BdG equations shows that there exists an intermediate 
critical coupling cJ-^ < oi\ where two real eigenvalues with 
m = coalesce and eventually become complex for a > a'j^ . 
The precise behavior of the complex eigenvalues with in- 
creasing a depends on the specific value of /3. A common 
feature is that complex eigenvalues as such again disap- 
pear before reaching the critical coupling a\. The latter 
is characterized by the fact that only one purely imag- 
inary eigenvalue persists for a > ai and signals a new 
bifurcation that leads to the appearance of axisymmetric 
vortex rings. 

As an example we consider the popular ratio /3 = 1/4 
for which we find ao — 10, c^^ = 32, a\ — 39. A static 
soliton exists for all a but is stable only for a < 10, 
while a static solitonic vortex appears as the most stable 
structure for a > 10. A new threshold occurs at a = 39 
where a static vortex ring emerges with energy interme- 
diate between that of a soliton and a solitonic vortex. 



in analogy with the situation on an infinitely cylindri- 
cal trap demonstrated in Fig. ^ Explicit illustrations of 
the currently calculated profiles are given in Figs. I8I9I for 
ct = 100 and /3 = 1/4 where all three types of static solitary 
waves are possible. For smaller (3 the critical couplings 
Q!o and ai increase indefinitely. But the corresponding ef- 
fective couplings 7o and 71, numerically calculated from 
Eq. (|B2|) applied for 2 = 0, converge to the asymptotic 
values 7o = 1.5 and 71 = 4 that are in agreement with 



N 




X y 

FIG. 9: Profile of a static solitonic vortex on a finite trap 
with a = 100 and P = 1/4, illustrated through density plots 
over two planes that are perpendicular to each other and both 
contain the symmetry (z) axis. Note that the solitonic vortex 
is not axisymmetric. 



the critical couplings obtained in the main text for an 
infinite cylindrical trap. It is also interesting to consider 
the special limit of a spherical {(3 = 1) trap in view of the 
experiment described in Ref . 0] ■ The critical coupling ao 
vanishes at /3 = 1 and thus the black soliton is never sta- 
ble on a spherical trap. On the other hand, the solitonic 
vortex now becomes degenerate with the ordinary vortex 
and is stable for all a > 0. The vortex-ring threshold is 
here predicted to occur at ai = 5. For the parameters 
of the actual experiment 7] we find a = 430 which is 
substantially larger than ai and eventually explains the 
observation of multiple vortex rings. 
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